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We investigate the statics and dynamics of spatial phase segregation process of a mixture of 
fermion atoms in a harmonic trap using the density functional theory. The kinetic energy of the 
fermion gas is written in terms of the density and its gradients. Several cases have been studied by 
neglecting the gradient terms (the Thomas-Fermi limit) which are then compared with the Monte- 
Carlo results using the full gradient corrected kinetic energy. A linear instability analysis has been 
performed using the random-phase approximation. Near the onset of instability, the fastest unstable 
mode for spinodal decomposition is found to occur at q = 0. However, in the strong coupling limit, 
many more modes with q ~ Kf decay with comparable time scales. 
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I. INTRODUCTION 



Recent realizations of twoEm and threeH component alkali Bose-Einstein condensates (-BEC's) in a trap provide us 
with new systems to explore the physics in otherwise unachievable parameter pregimes.uffl. Dramatic results have 
recently been observed in the phase segregation dynamics of mixtures of RlxrEI and NaB gases. Periodic spatial 
structures were found at intermediate times which then recombine at a later time. 

Phase segregation phenomena have been much studied in materials science and these can be-.undcrstood using clas- 
sical mechanics. Spatial modulations have also been-observed, for example, in AINiCo alloys.LJ These were explained 
• i-h . in terms of a concept called spinodal decomposition.!! When a system is quenched from the homogeneous phase into a 
broken-symmetry phase, the ordered phase does not order instantaneously. Instead, different length scales set in as the 
domains form and grow with time. For the BEC's, however, quantum mechanics play an important role. It has been 
O ; shownQ that it is possible to have an analogous spinodal decomposition, which manifests some of the phenomenology 
'— 'j including a periodic spatial structure at an intermediate time that is now determined by quantum mechanics- The 
time scale provides for a self-consistent check of the theory and is consistent with the experimental results.13 The 
growth of domains at later times, is now determined by quantum tunneling and not by classical diffusion. 
00 . Recently, it became possibleE3 to cool a single component system of about a million 40 K fermionic atoms in a 
£SJ magnetic trap below the Fermi temperature, Tp, leading to the realization of a spin-polarized fermion gas of atoms. 
Similar to electrons in a solid, the dilute gas of atoms fills all the lowest energy states below the Fermi energy, Ep. 
The transition to this quantum degenerate state is gradual as compared to the abrupt phase transition into a Bose 
condensate. For single component fermionic systems, however, the equilibrium is difficult to achieve aSr-the s-wave 
elastic collisions are prohibited due to Pauli exclusion principle. In the experiments of DeMarco and JinE3, this was 
circumvented by using a mixture of two nuclear spin states of 40 K atoms for which s-wave collisions are allowed. 
One of the manifestations of quantum mechanics was the nature of momentum distribution which differed from the 
. j—{ ' well known classical gaussian distribution. This system corresponds to the weak coupling limit in which the-physical 
properties are close to those of a non-interacting fermion gas. The other system which is being exploredHl is the 
gas of 6 Li atoms. Mixtures of fermions interating with the Coulomb interaction have been studied in the context 
of the electron-hole fluidaI3. For fermions mixtures on a lattice site interacting with the Hubbard Hamiltonian, the 
■ partial phase segregation leads to antiferromagnetism. Thermodynamic properties as well as density and, .momentum 
distributionSj-of spin-polarized fermionic gas of atoms in a harmonic trap have been studied in recent yearsoli5E£l. Butts 
and Rokhsarli3 have obtained universal forms of the spatial and momentum distributions for a single component spitu 
polarized non-interacting fermion gas using the Thomas-Fermi (TF) approximation, whereas Schneider and Wallisc£l 
have studied the effects of shell closure for small number of atoms, similar to the nuclear shell model. Bruun and 
Burnetto have studied an interacting fermion gas of 6 Li atoms .which have a large negative scattering length. Such an 
interaction could also lead to the possibility of superfluid stateliS in these systems. In the present paper, we consider 
mixtures of these new finite systems of ultracold fermionic atoms with a positive scattering length in the limit of both 
weak and strong coupling and explore the equilibrium and non-equilibrium quantum statistical physics using the TF 
approximation, Monte Carlo simulations, and the random phase approximation. 

In section II we present the equilibrium static properties of mixtures of fermionic atoms in different parameters 
regimes using both the TF and the Monte Carlo simulations. In section III, we study the dynamics of phase segregation 
of such mixtures using a linear stability analysis. Finally, conclusions will be presented in section IV. 
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II. STATICS 



We first start with the statics of a two component fermion gas of atoms with masses m\ and iri2 and particle 
numbers N\ and N 2 . This is assumed to be confined in an azimuthally symmetric harmonic trap with radial and axial 
frequencies u> and Xlo, respectively which are considered to be the same for both the components. Unlike the electron 
gas in matter, the fermion gas of atoms is neutral and dilute. Therefore, significant interactions between atoms are 
only short-ranged and that would be responsible for any phase segregation in the system. In the long wavelength 
limit, the system can be well described by the density functional theory and the total energy can be written as 



E = J £ EoM + ,gpi(r)p 2 (r)]dr. 



(1) 



Here Eo a = j^-T a (r) + ^m lT uj 2 (x 2 + y 2 + A 2 z 2 )p CT (r) is the non- interacting part of the energy density and p CT (r) is 
the particle density of the component a = 1, 2 with J p&(r)dr — No-. The interaction term has been approximated by 
the contact potential gS(r — r'). g is related to the scattering length a by g = 2irh 2 a/m : with m = mimzKmx + m 2 ). 
In accordance with the experiments, we take a to be positive and consider only the s-wave scattering. Therefore, the 
contribution to the interaction term is non-zero only when the species are different or are in different hyperfine states 
as in experiments. From the Pauli exclusion principle, there is no contact interaction between particles of the same 
species (spin). In a more general treatment including p-wave scattering there would be additional terms involving 
interaction between identical species also. But these are small, and thus neglected. 

For the kinetic energy density r CT we use a local approximation including the first and second derivatives of the 
particle density, 



. W = |(6^ W V» + ^^W + lvV(r). 



(2) 



The first term represents the Thomas- Fermi (TF) approximation to the kinetic energy. The second term is ^ | V v / /v | 2 
and represents the gradient correction to the kinetic energy. The integral of the third term extended to infinity vanishes, 
and thus it will not be included in the calculations. The Monte-Carlo results confirm that the gradient term is at 
least 2 orders of magnitude smaller than the TF term, but this term is important in that it can break the symmetry 
of the ground state and lead to asymmetric states with a lower energy . 

Without the interacticp. term in (1), the system behaves in the same fashion as the one component system for 
which Butts and RokhsarLJ obtained Ep to be related to the total particle number N by Ep = hui(6XN) 1 ^ 3 . Defining 
Rp = (2Ef /muj 2 ) 1 / 2 (giving the characteristic size of the gas), and Kp — (2mEp /ft 2 ) 1 / 2 (momentum of a free 
particle of energy Ef), they calculated the density profile at T=0 to be given by 



Pn 



-interacting (r) = PO [l ~ r 2 / R%] 
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(3) 



with r = x + y + A z , po = 8NX/n Rp = K f /6tt . In the TF approximation, the trapping potential can be 
treated to be locally constant and we can define a local Fermi wavevector, kp(r) so that Ep = h k%(r) /2m + V(r), 
and the density at T — can also be written as p n on-intoracting(r) = kp(r)/6ir 2 . 

We now examine the properties of the mixed (two-component) interacting system and will show how the repulsive 
interaction modifies this non-interacting density profile as well as other properties of the system. The strength of the 
coupling, which controls the phase segregation, depends on the dimensionless parameter which is the ratio between 
the interaction and the kinetic energies, namely ffPiP2/[-7g-(67r 2 ) 2 / 3 (p 1 ^ 3 jm\ + p^ 3 / 77 ^)]- In the simple case of equal 
masses (mi = m 2 = m) and densities (pi = p 2 = p) of the two components, this simply scales as aKp. This means 
that the coupling would be stronger if a or the density is large. Also as Ep is proportional to the frequency of the trap 
at constant N (a higher frequency leads to a larger separation between the levels), the coupling would be large for 
higher frequencies. From now on, to measure the strength of the interaction, we will use the dimensionless parameters 
Ccr = Kp a a/n, where Kp a = (2m cr p cr ) 1 / 2 jfi , or in the case of equal chemical potentials, just c = Kp a/ir. 

For a general two-component system with chemical potentials p,\ and ,u 2 , the ground state is obtained by minimizing 
the thermodynamic potential fl — E — J(pipi + ,u 2 p 2 )dr. This leads to the following system of equations: 
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Similar to the one-component case, one can rewrite the above in a dimensionless form by introducing for each of 
the species a, the following quantities: R a = [2// CT /m CT o; ]s , p a Q = Kp a /6ir 2 , Q a = gpsol '/V, an d n a (r) = Pa(r)/p a Q. 
Here a = 3-cr. If one neglects the smaller terms containing derivatives of p (the TF limit), one obtains the following 
algebraic equations satisfied by the dimensionless densities n\ and n 2 for any coupling strength Q a : 

n\ /3 = 1 - f 2 /R\ ~ Q x n 2 

n 2 /3 = l-f 2 /R 2 2 -g 2ni . (6) 

We see that the effect of the additional Q„n s term, i.e. the interaction, is to deplete the regions where n s is highest 
(without necessarily leading to a phase segregation). 

When there-is phase segregation, the interface energy is proportional to the square root of the coefficient of the 
gradient termEZI and it often serves to distinguish different configurations. In that case, their effect cannot be neglected 
and these are included in the Monte Carlo simulations. We next discuss some special cases in the TF limit. 

A. TF limit: Similar densities: (pi — ^2) for any coupling 

To simplify the notations, we will use: p,\ = fi 2 = A 4 ; Ri = R2 = R) Q\ = Qi = Q ■ In this case, three solutions 
to Eq. (6) will correspond to n\ — n 2 , of which only one is physical with ri\ > 0. If a solution n 2 = f(n\) exists, 
by symmetry, the other one is necessarily n\ = f(n 2 ). These solutions with n\ 7^ n 2 can be obtained numerically. 
The real solutions are plotted in Fig. [I], where the ri\ = n 2 solution is referred to as "Sym" , and the other conjugate 
(asymmetric) solutions are referred to as "Al" and "A2" . Below we discuss these solutions in the weak and strong 
coupling limits. 

1. Weak or intermediate coupling regime 

In this case we look for symmetric solutions (rii — n 2 = n). Equation (g) then reduces to (dropping the subscripts): 

n(r) 2 / 3 = l-f 2 /R 2 -Gn(r), (7) 

which can be solved easily numerically to give the density profile of the non-segregated phase. It is possible to show 
that after proper rescaling, the result for all coupling strengths and at any point can be summarized in a single universal 
curve in Fig. [j]. If n(r) is a solution to Eq. (|7|), then Af = nQ z versus V = [l — f 2 /i? 2 ] Q 2 is the universal function 
of Fig. satisfying M 2/3 + N - V = 0. For small couplings and near the boundary (V 0; Af 2/3 > Af = V 3/2 

), this curve is a power law and in fact tends to the non-interacting density n(r) ra [l — (x 2 + y 2 + A 2 z 2 )/i? 2 ] 3 ^ 2 . 

2. Strong coupling regime 

The above situation, however, can not be always sustained. In the strong coupling limit, we can have phase 
segregation (n\ 7^ n 2 ), and one needs to go back to Eq. (^) which now admits lower energy solutions that are not 
"permutation symmetric" : 

nl /3 + Gn 2 = l- (x 2 +y 2 + \ 2 z 2 )/R 2 & < /3 + Af 2 = V N 2 = (V - M 2 ) 3 

n 2 2 /3 + g ni = 1- (x 2 +y 2 + \ 2 z 2 )/R 2 N^' 3 + M = V «*■ Af = (V ~M) 3 i (8) 

where we used the same simplifying notations as before. As previously mentioned, the symmetric solution Af± = Af 2 
always exists. This can be exploited to reduce the above equations to a quadratic equation, which is analytically more 
transparent. 

Subtracting the above equations from each other and dividing out by Af± — Af 2 , we obtain, 

A - ! + Af 2 - {V - Af 2 ) 2 + (V -Afi) 2 + {V - Af 2 )(V - Af x ). (9) 

This quadratic equation can be solved for Ai in terms of Af 2 . 

The solutions will all be axially symmetric in that they are functions of f 2 only. In actuality, the axial symmetry 
can also be broken, but we do not find it here since we neglected the terms in gradient of the particle density in 
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the kinetic energy. The broken symmetry solutions will be discussed in the subsection E where we present results 
obtained from the Monte Carlo simulations incorporating these terms. In Fig. |], the solutions with n\ 7^ ri2 can be 
seen in the limit of small reduced distance and large V . The bifurcation point where these solutions start to occur, 
corresponds, from numerical results, to V c ~ 0.741, and Af c — nQ 3 0.296, which separate the strong coupling regime 
from the weak one. In both figures, the symmetric solution is drawn with solid line, and the asymmetric ones with 
dashed lines. Actually, at the bifurcation point, we have exactly Qnz = 2/3 as will be shown in the TF linear stability 
analysis section below. Since G 2 =V/(1— f 2 /R 2 ) > V, the smallest coupling Q c for the unequal solutions to occur 
satisfies Q c — \fW c . Since Q = (4/3)Kpa/7r, we find a critical dimensionless coupling c = (Apa/7r) c ~ 0.646. We shall 
come back and compare this value with that obtained with a different approach. 



B. TF limit: Very different densities: (pi S> P2) for any coupling 



One can also treat the case where one of the species is a minority (pi ^> P2). If we assume p\ = \ 2 p2 1 then 
Ri = Ai?2; Kfi = \Kf2]Pw — A 3 p2o; Gi — A Gi, and n a ~ 1. The density distribution of the majority species will 
be weakly perturbed. Referring to Eqs. (Q), one can see that the coupling Gi — gp2o/Pi becomes very small and 
maybe neglected. Thus a good approximation is to assume p\ ~ jO non _i n t e racting- The G2 term in the second equation, 
however, is a large quantity, and will strongly affect the particle density U2- Therefore, 
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(10) 



In the presence of the majority species, the number of atoms of minority species will be much less than their non- 
interacting counterparts with the same chemical potential. As we can see from the above equation, their number, 
even at the origin is reduced by a factor of (1 — G2) 1 ■ We find that for a large enough G2 the density A/2 is depleted 
from the center (see also Fig. [j]b, curve A2). 



C. TF limit: linear instability analysis 



We next study the fluctuations of the system about its equilibrium configuration in the TF limit by expanding 
the thermodynamic potential Q upto second order in the particle density variation bp about its minimum which was 
computed above. The sign of the second derivative of ft will decide the stability of the symmetric phase. A phase 
segregation occurs when the Hessian (second derivative matrix) ceases to be positive definite. If the transition is first 
order, it would have already occurred before reaching a negative second derivative. The second derivative from Eqs. 
(3) and (4) is just a 2 x 2 matrix: 

d 2 n h 2 2, 9n 2 _i r r , 

-o( 67r Y p° <W +s(i-<W)- (11) 



dpadpa' 1rn a 3 

The phase instability criterion thus becomes uj- = where w_ is the smallest eigenvalue of the Hessian matrix; 
implying: 

H( 67r 2)§ ( pi p 2 )-§ = g&JLl n -l/a = gi £( pi=pa ) (12) 



2 v / mim2 3 po 3 

Thus, in the symmetric case (/ii = P2] Pi = Pi), the instability will first occur locally at the point where the relation 
Af 1 ' 3 = Gfis — 2/3 is satisfied. This implies that J\f — 0.296, which is exactly the critical A/" c obtained earlier from a 
different analysis. These two instabilities occuring at the same point suggest that, within the adopted model (TF), 
the transition might be of second order. 



D. Possibility of density modulation instability 

Similar to the electron gas which has several kinds of instabilities such as ferromagnetism, antiferromagnetism, 
charge density wave, superconductivity, etc... these two-component systems might also exhibit other types of insta- 
bilities. To investigate them, we will assume the homogeneous case (w = 0) as the analysis can be made simpler by 
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using the Fourier decomposition of the density. To get some understanding of the nonuniform systems (such as in a 
trap), one can assume in a semiclassical approximation, that the Fermi momentum depends on the position, as before. 

The density for the species a can be written as the sum of its Fourier components: p a {v) = p a + X) q ^o P°i eJqr > 
with p a 3> Pa-q- Substituting this expression in the thermodynamic potential f2, expanding up to second powers of 
p a q, and minimizing with respect to the Fourier components, we obtain: 



dp a q 
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(14) 



Assuming 6ir 2 p~a = k^ (note that in the presence of interactions, the average density and Fermi momentum, which 
we denote here by Pa and k a respectively, are different from their non-interacting values), the above equations are 
simplified to: 
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(15) 
(16) 



It is clear from the above equations that if a = then p a q = is a solution (uniform density if no coupling) . For 
a > 0, we have p a q and p S q of opposite signs for all q. This means that there is phase segregation for repulsive 
couplings. Furthermore, if a < 0, there will be density modulation in the small q limit (the functional we considered 
is valid in the long wavelength limit). We shall return to this point in section III where the dynamics are treated. 

One can also note that the transition points of Eq. ^ and previously studied Eq. |^ are the same (in the u> = and 
pi = pi case), since they are derived from the same functional. Indeed from the positive-defmiteness of the functional 
Q in this representation, one obtains that the transition occurs for kaa/n = 1/2. Inserting this critical value into 
Eq. [IB], one finds the relati on b etween the non-interacting Fermi wavevector Kf and the interacting one k a at the 
transition point: Kp = k a --< ] /5/3 which then implies 



= ^/^ = ^1-0.645, 



(17) 



which is exactly the same value as obtained from the numerical result of the previous section. 



E. General case: Monte Carlo results 



The density distribution that extremizes the energy functional in Eq. (1) can be obtained by a Monte Carlo 
simulation with a weighting factor exp(—E/T) for a parameter T that is sufficiently, low. This is basically the 
simulated annealing method and has been exploited successfully in earlier treatmenttj of the corresponding Bose 
system described by a Gross-Pitaevski functional. 

We approximate the volume integral of the energy functional by a discrete sum. Using the scaled radius f, we 
sample a lattice inside a sphere of diameter 2R consisting of 40 sites along the diameter, making a total of 33398 sites. 
The derivative term is approximated by a finite difference. For simplicity, we show here results for the case when the 
two components have the same mass. 

We first show in Fig. [2] the density profile of component 1 as a function of x and y for z=0 for the weak coupling case 
with no phase segregation. The values of different parameters were chosen to be w = 135 x 27rrad/sec, a = 135aBohr, 
A = 0.14, and -/V±i = ^2 = 10 6 (pi = P2 = 1-626 x 10~ 29 J); roughly corresponding to the experimental parameters of 
the 40 K systemEj. In these experiments, we estimate c = Kpa/n = 0.032, Rp = 26 pm, and Q = 0.042. The density 
profile for component 2 is the same and hence is not shown. 

In the limit of strong interaction, phase segregation starts, and as mentioned earlier, the system can now also break 
cylindrical symmetry. This happens when Kpa is large enough, which in turn can be achieved with only large Kp, 
only large a, or both. To illustrate this, we show in Fig. Stlie density profiles for components 1 and 2 for the case of 
only large a with a = 30000 aeohr, Pi = P2 — 1-86 x 10 J, and u> = 300rad/sec. In this case, c = Kpa/n = 2.39, 
Rp = 25pm, and Q — 3.19. 



5 



For the case of both large Kp and a, we show in Fig. |J the density profiles for a — 3000aB O hr, Mi = M2 = 
2.762 x 10~ 29 J, and w = 2000rad/sec. This corresponds to c = 0.92, R F = 14.3 /an, and = 1.23. The difference 
in the densities of the two components shows that the largest change occurs near the center where the density is 
maximum. 

It is to be further noted that for this case, the density distribution is still quite cylindrical but there is a slight 
asymmetry, as we can see from the graph of the difference. This asymmetry becomes more pronounced as the 
interaction is increased further. In Fig. ^ we have shown the results of simulations with larger a. The density 
profiles were calculated for a — 4160aB O hr, Mi = M2 = 1-626 x 10 _29 J, ui — 6000rad/sec. This corresponds to 
c = 0.98, R F = 3.67 /xm, and Q = 1.31. 

As discussed earlier, phase separation can also occur when Ni >> N2. As an illustration, we show in Fig. || the 
density profiles for components 1 and 2 for the case a = 104aBohn Mi = 2.6016 x 10~ 26 J, p 2 = 4.336 x 10~ 26 J and 
u) = 1600000 rad/sec. 

The density of component 2 is small and therefore, its noise is also substantially higher. One can clearly see the 
density depletion of component 2 at the center. 



III. DYNAMICS 



We next turn our attention to the issue of dynamics. For the classical and boson spinodal decompositions, the 
fastest unstable mode occurs at a finite wave vector. We ask if a similar situation occurs for the fermion case. We 
found that the fastest unstable mode occurs at wavevector q = at the onset of instability. For stronger coupling, 
many modes with q ~ Kf decay with comparable time scales. We now describe the details of this linear stability 
analysis. 

The energy functional (Eq. (1)) which was approximated with a local kinetic energy depending on the density and 
its derivatives is only good in the long wavelength limit. Due to this approximation, we found that the instability has 
a local character and occurs first in regions of high density. Here we will perform a linear instability analysis in the 
random phase approximation (RPA) to improve upon this local picture. The linear susceptibility \ is defined as the 
response of the particle density to an external potential V which could also be tr-dependent: 



5p a {v) = V / dr' x<7a ,(r,r')V^(r'). (18) 



E 

er'=l,2 

Here V tot is the total self-consistent field and is the sum of the external field and that due to, the interaction: 
V^ ot = V a + gSpa- The bare response x<ya can be obtained from the usual Lindhard expression^. Since there is 
no term in the Hamiltonian that interchanges the species 1 and 2, off-diagonal terms of the susceptibility are zero 
(Xi2 = X21 — 0). Taking the above into consideration, Eq. (|l8j) can be written in the following matrix form: 
5p = x(V + G5p), leading to 5p = [1 — xG]~ 1 x V, where the 2x2 matrix G has as its diagonal elements and g as 
its off-diagonal elements, and x is diagonal. Consequently, an instability will occur when the following determinant 
becomes zero: 

Dct|l- X G| = l-. 9 2 xiiX22 = 0. (19) 
In the case where the densities are equal, xn = X22 = X, the two eigenmodes are calculated as: 

Spi+6p 2 = (i- X gy 1 x(v 1 + v 2 ) (20) 
5 P i - S P2 = (1 + xgr'xiVi - V 2 ). (21) 

The first mode corresponds to a density fluctuation, and the second mode Spi — Sp 2 represents the phase separation 
instability in which we are interested. The response corresponding to this mode is given by e(q,w) = [1 + gx(q,w)]. 
The instability decay time v~ x is determined from the formula e(q, iv) — 0, since, in this case, any infinitesimal 
external potential will lead to a large change in the density. There exists a q — qo such that v(qo) is largest. This 
determines the spinodal wavevector of the fermionic system as it indicates the mode with fastest growth. In what 
follows, we will be treating the constant external potential problem where the Fermi momentum is k. For the confined 
case, one can consider k to be a local function related to the density by k(r) — [67r 2 p(r)] 1 / 3 . From the Lindhard 
expressionES for x (real frequencies), we obtain, after correcting for a spin degeneracy factor of 2, the corresponding 
dimcnsionless response x = — 47r 2 ?i, x(<?> *^)/wfc for imaginary frequencies: 

X (q^) = l + -(l + (v/q) -( g /2))Log[ (1 _ g/2)a + (t//g)a ] (22) 
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-iK 1 ^ 1 ^- 1 '^ 1 )- 

Here q is in units of k and v, in units of K/2E — m/Tik 2 . The three-dimensional plot of x as a function of q and v is 
shown in Fig. |?]. 

The equation e(q, w) = [1 ± gx(q, w)] = implies that the instability points for phase segregation with a repulsive 
interaction [g > 0) and that of density modulation with an attractive interaction are exactly the same within RPA. 
This is also in agreement with the analysis of section II E where it was shown that "magnetic" instability occurs for 
repulsive interactions, and "density wave" instability may occur for attractive interactions. Although the susceptibility 
can be both negative or positive, for a coupling of fixed sign, one should only consider the physically correct situation. 
In our case, for positive g, only the "magnetic" instability, i.e. \ — ~~ Vs should be considered. 

Now since gx = —xka/it, the instability condition implies c\ = 1 where c = ka/n. The maximum of x is obtained 
for q — > and lo — > where it tends to 2. From this result, we arrive at the conclusion that there is no solution to 
e(q, iv) — for c < 0.5 and no instability develops. For larger values of c, the plane z = 1/c intersects the surface of 
X on a curve which is displayed in Fig. |8|. The inverse decay time lo as a function of the wavevector in units of k is 
shown in this figure. As can be seen, the fastest unstable mode occurs at wavevector q = and lo = at the onset 
of the instability (c = 0.5) in agreement with Eq. [I?] previously derived. Indeed the instability calculation derived in 
the previous section focused on the long wavelength aspect of the problem. 

For stronger couplings, many modes with q sa k decay with comparable time scales of the order of h/Ep, but those 
with shortest timescales (i.e. largest lo) prevail. |— ,|— , 

In the really strong interaction limit, further phase separation can take place either via tunnellingtaEJ or via 
quantum motion of the domain walls. We hope to investigate this further in the future. 

The behavior of the wavevector of instability is similar to that of the classical spinodal decomposition, which we 
briefly recapitulate here. The current J can be related to the free energy F by Fick's law: J = cVF for some constant 
c. After the onset of instability, F = (—A + Bq 2 )8p q . As one goes from the onset of instability, A starts to become 
non-zero. In addition, there is the particle conservation equation — dtp = V • J. Combining the above two equations, 
we obtain ito5p q = cq 2 (—A + Bq 2 )8p q . The fastest mode occurs at a wavevector q c = yj A/2B. Thus at the onset of 
instability, q c = 0. q c becomes larger as one goes away from the instability point. 

IV. CONCLUSION 

In conclusion we have investigated the statics and dynamics of the spatial phase segregation process of a mixture 
of fermion atoms in a harmonic trap using the density functional theory and the random phase approximation. As 
the coupling starts to increase, even with the same chemical potential, equilibrium distribution with unequal densities 
starts to appear, which quite often do not exhibit axially symmetric correlations. Similar to the classical and Bose 
spinodal decomposition cases, the fastest mode for the initial phase segregation occurs at a finite wave-vector. The 
condition of instability corresponds to a large interaction, which may be achieved experimentally with the atoms close 
to a Feshbach resonance. 

The instability calculation for the phase segregation phenomena discussed here is related to the instability calculation 
for the antiferromagnetic transition of the electron gas. In the electron gas, this is enhanced when there is nesting of 
Fermi surface such as in Cr or in one dimensional materials. The transition always stops after the 2Kp instability 
due to the long range nature of the Coulomb interaction, and no further "segregation" takes place. 

An interesting situation is the one dimensional trap as it would exhibit a much stronger instability. In mean 
field, the one dimensional density difference response function e(2Kp) = 1/[1 + Kpa Log(T/Ep)] is logarithmically 
divergent at zero temperature. The transition temperature occura.|a± T c = Epe~ 1 / KF ' a . One dimensional trap, which 
can be realized for small values of A, has been extensively studiednEa and we expect a higher tendency towards phase 
segregation in that case as well. 
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FIG. 1. Top: Dimensionless density versus dimensionless radius f/R for Q — 1. One of the asymmetric solutions (A2) is 
depleted at the center while the other one has a large concentration. For f/R larger than 0.51 both asymmetric solutions join 
the symmetric density profile. The sharp features around this point are due to the neglect of the gradient terms. Bottom: 
Universal curve of rescaled density Af — nQ 3 versus rescaled distance from the border V = (1 — f 2 /R 2 ) Q 2 , valid for all coupling 
strengths Q. Note that < V < 1, and for the symmetric case A/max = 0.43 (f = or V — 1). 

FIG. 2. Snap shot of the density profile at z=0 as a function of x and y in the weak coupling limit c = 0.032. 

FIG. 3. Snap shot of the density profile of components 1 (top) and 2 (bottom) at z=0 as a function of x and y in the strong 
coupling limit c = 2.39, u> = 300 rad/sec. 

FIG. 4. Snap shot of the density profile of components 1 and 2 and their difference at z=0 as a function of x and y in the 
strong coupling limit (c = 0.92, lo = 2000 rad/sec). 

FIG. 5. Snap shot of the density profiles of components 1 and 2 at z=0 as a function of x and y in the strong coupling limit 
(c = 0.98, w = 6000 rad/sec). 

FIG. 6. Snap shots of the density profiles at z=0 as a function of x and y for ci = 0.98, C2 = 1.27, u> = 1600000 rad/sec. 
Density 2 is depleted in the central region. 

FIG. 7. Surface plot of the positive part of the reduced Lindhard susceptibility (x) as a function of q/k and the imaginary 
frequency. 
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Fig. 2, K. Esfarjani et al. 
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Fig. 3 (top:density 1), K. Esfarjani et al. 



12 




Fig. 3 (bottom:density 2), K. Esfarjani et al. 
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Fig. 4 (middlerdensity 2), K. Esfarjani et al. 
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Fig. 4 (bottom:density difference), K. Esfarjani et al. 
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Fig. 5 (bottom:density 2), K. Esfarjani et al. 
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Fig. 6 (top:density 1), K. Esfarjani et al. 
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Fig. 6 (bottom:density 2), K. Esfarjani et al. 



20 




21 



1.4 
1.2 



0.8 
0.6 
0.4 
0.2 




Fig. 8 (contours of chi bar), K. Esfarjani et al. 
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